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Abstract 


The first part of the present thesis describes an experimental technique for the inverse deter- 
mination of thermal conductivity in a 2D square domain based on a finite-difference based 
numerical algorithm, which uses transient temperature data as the input.The algorithm is 
tested for a known conductivity material such as mild steel.lt is found that the algorithm 
is highly sensitive to the measurement errors in the input data, a fact not unexpected in in- 
verse problems which are essentially ill-posed in nature. In addition to the above mentioned 
experimental work, a detailed stability analysis is also performed in order to find the maxi- 
mum permissible measurement en*or for such problems using simulated temperature profiles 
containing random eiTors. Interestingly, the study reveals that the maximum permissible 
measurement error is high at the early time and decreases exponentially with time. This 
means that early time temperature profiles should be used as input data. However, even at 
early times maximum allowable eiror in the temperature measurement is quite small which 
explains why. the present experimental input data did not give correct conductivity value for 
mild steel.Future effort is needed in the direction of improving the expenmental set-up so 
that high precision temperature measurement is possible. 
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Nomenclature 


Cp specific heat of the material at constant pressure 
h convective heat transfer coefficient 
k thermal conductivity of the material 
kmax Maximum conductivity value in the matrix 
kmin Minimum conductivity value in the mati’ix (”]f ) 
kth Theoretical conductivity value 

kouput Conductivity value matrix obtained from inverse problem 

kinput Conductivity value matrix given as input to get simulated temperature data 

L length (rn) 

n number of grid points in x direction and in y direction 
t time (second) 

T temperature [°C) 

Tjn Mean temperature value {'’C) 

Tom Mean temperature at 0 th time (‘’C) 

T' Maximum noise that can be allowed theoretically (°C) 

Maximum noise at 0 th time that can be allowed theoretically (°C) 

To Temperature at 0 th time (Tom. + 

Zl, Z2, Z3, Z4.l,ZA.2 Stability parameters 
6 Error introduced in simulated temperatures {‘'C') 

5max Maximum allowable error in simulated temperatures to get reasonable value of 'k’ by 
using inverse simulation (fC) 

So Maximum allowable error at 0 th time {°C} 
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Greek Symbols 


a relaxation factor 

e small value of temperature measurement error 
p density of the matenal (fcp/m^) 

Subscripts 

1 at first of the two times at which input temperature is taken, also top or left boundary. 

2 at last of the two times at which input temperature is taken, also bottom or right boundary. 
i space index in x direction, also at t=0 

j space index in y direction, 
oo ambient 


Superscripts 

p present time 
p + 1 future time 
' first derivative 

" second derivative, also per unit area 


Special Symbol 


A increment 



Chapter 1 


Introduction and Literature Review 


1.1 Introduction 

The use of inverse heat conduction techniques for the determination of thermal prop- 
erties such as thermal conductivity and heat capacity of solids or the estimation of surface 
condition such as temperature and heat flux by utilizing the transient temperature measure- 
ments taken within the medium, has numerous practical applications. 

In the past decade, rapid advancement has been made in the field of materials research. 
This is primarily due to the demand for advanced material witli low weight or the ability 
to withstand high temperatures required in the nuclear, electronics and aerospace industnes. 
Precise knowledge of thermophysical properties for these materials is essential in many ther- 
mal systems. Specifically, an accurate prediction of thermal conductivity is imperative to 
achieve thermal control system. 

The direct measurement of heat flux at the surface of a wall subjected to a fire, at 
the outer surface of a re-entry vehicle or at the inside surface of a combustion chamber 
is extremely difficult. In such situations, the inverse method of analysis, using transient 
temperature measurements taken within the medium can be applied for the estimation of 
such quantities. 

An inverse problem is much more difficult to solve than the direct one. The reason for 
this is that it is usually ill-posed, i.e., it is very sensitive to measurement errors. An excellent 
discussion of difficulties encountered in inverse analysis is well documented in Beck et al 
(1985). 

The present work involves the estimation of thermal conductivity by an inverse tech- 
nique using finite-difference. Therefore, the literature review that follows concerns only the 
papers and texts associated with thermal conductivity measurements. 



1,2 Literature Review 


To date,vai'ious methods have been developed for analysis of the inverse heat conduc- 
tion problem involving the estimation of thermal conductivity from measured temperatures 
inside the materialBeck and Arnold(1977) determined the thermal conductivity by mini- 
mizing the eiTors between the measured and calculated temperatures in a least squares sense 
Alifanov and Mikhailov (1978) used the conjugate gradient method to solve the non-linear 
inverse thermal conductivity problem.Chen and Lin (1981) developed a pulse-spectrum tech- 
nique (PST) based iterative numerical algorithm for remote sensing of the thermal conduc- 
tivity of a non-homogencous material for the onc-dhncnsionaJ (ID) casc.it is found that PS f 
does give excellent results and is more robust in solving the inverse problem of the diffusion 
equation than that of wave equation. Jarny et ai. (1986) employed an Output Least Square 
Method (OLSM) to determine thermal conductivities of materials. 

Tervola (1989) discussed the use of finite element method in conjunction with Davidson 
Fletcher-Powell method for determination of thermal conductivity of a homogeneous mate- 
nal from the measured temperatures at certain ( finitely many ) points.The work concerns the 
situation where the thermal conductivity is dependent on temperature. Jaimy et al, (1991) em- 
ployed a general optimization method using adjoint equation for solving multi-dimensional 
inverse heat conduction problems. Chen and co-workers (1996) employed the hybrid scheme 
of the Laplace Transform technique and the central difference approximation to estimate 
temperature dependent thermal conductivity from temperature measurements inside mate- 
rial at an arbitrary time.Lin and Chang (1997) proposed a numerical algorithm to estimate 
the tliermal conductivity of a homogeneous material from boundary temperature measure- 
ments.Dowding et al. (1996) presented a laboratory method to measure the thermal prop- 
erties of a carbon-carbon composite material that is characterized by an orthotropic thermal 
conductivity and isotropic volumetric heat capacity. Alencar Jr. et al. (1998) presented a gen- 
eral solution for two-dimensional boundary inverse heat conduction problem by using the 
conjugate gradient method of minimization together with an elliptic scheme of numerical 
grid generation, Simulated measurements are used to illustrate the application of the present 
approach to the solution of an inverse problem of practical interest. 

The literature on the inverse determination of thermal conductivity from measured 
temperature data by finite-difference method are only a few. Lam and Yeung (1995) em- 
ployed a first order finite-difference method to determine thermal conductivity in a one- 
dimensional (ID) heat conduction domain using Cartesian coordinates. In a subsequent 
paper, Yeung and Lam (1996) have successfully extended their earlier results,examining the 
feasibility of using a second-order finite-difference technique to determine the thermal con- 
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ductmty in a one dimensional (ID) heat conduction domain Recently, Ghoshdastidar and 
Ray (1998) presented a second-order finite-difference procedure for the inverse detennina- 
tion of thermal conductivity in two-dimensional (2D) square and cylindrical domains using 
available transient temperature data at discrete grid points. The procedure is capable of pre- 
dicting constant, spatial and temperature-dependent thennal conductivities. For 2D square 
and cylindrical domains, the estimated thermal conductivity is verified for a constant con- 
ductivity material such as steel using numerically simulated transient temperature data.For 
2D square domains, the results have also been tested against five bench mark solutions for 
ID domain in which thermal conductivities are either constant or linear/non-linear function 
of space and temperature by forcing a 2D domain to behave like a ID domain.The close 
agreement between the current results and the exact solutions confirm accuracy and effec- 
tiveness of the proposed finite-difference technique. It may be noted that the advantage of 
using finite-difference method is that the conductivity function can be obtained by solving 
the system of linear equations arising out of the discretization of the governing partial dif- 
ferential equation. Therefore, no prior information is required on the functional form of the 
thermal conductivity. 


1.3 Objectives of the present study 

The objectives of the present study are as follows. 

1 . To extend the work of Ghoshdastidar and Ray (1998) to the special case of a 2D square 
domain of a known constant conductivity material (mild steel) in which no heat is gen- 
erated.Thus, the transient temperature data are created by obtaining numerically sim- 
ulated temperature for the situation when the body is heated to a uniform temperatuie 
(Ti) and then is allowed to cool by losing heat to the ambient. 

2. To develop an experimental set-up which will produce transient temperature data that 
will be used as the input to obtain the thermal conductivity of mild steel in a 2D square 
domain. 

3. Since, inverse heat conduction problems are very sensitive to measurement errors, an 
effort is made to quantify a realistic error in conductivity by adding random errors to 
the simulated temperature inputs 

4. To calculate conductivity using filtered Tactual where Tactual = Tea^^ct+measurement 
error 

5. To do a stability analysis based on 3 and 4. 
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1.4 Organization of the Thesis 


The present thesis has been organized in the following manner. 

> In Chapter 2. the finite-difference technique of Ghoshdastidar and Ray (1998) is dis- 
cussed with respect to the present problem. 

t> Chapter 3, details the experimental set-up used to obtain transient temperature data 
required for the numerical algorithm and gives estimated conductivity values of mild 
steel using actual and filtered temperatures. The sources of experimental errors and 
possible ways to alleviate them are also mentioned. 

o Chapter 4 presents the results of the stability analysis based on numerically simulated 
temperature inputs (Tactual) with the addition of the random measurement errors as 
well as that based on Tactual after filtenng using Gram Orthogonal polynomial method 
described in Al-Khahdy (1998). 

> Finally in Chapter 5, the conclusion drawn from the study and the scope for the future 
work are indicated. 
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Chapter 2 


Problem Formulation 


2,1 Problem Statement 

Consider the problem’ of unsteady heat conduction in a square plate (L x L) which 
is initially at a tempeiature, T^. Suddenly, all four sides are exposed to a cool ambient 
at Too. The problem can be assumed to be a two-dimensional one as the plate is thin so 
that the thermal picture in the plane normal to the Z-direction is more or less identical. The 
physical problem and computational domain is pictorially depicted in Fig.2. 1 . The simulated 
temperature profiles at various times have been computed using Explicit Finite-Difference 
scheme ( Appendix ’A’ ). 


2.2 Governing Differential Equation 


The governing differential equation for this case is, 

~ u{Q^^§^\jL.9kdT,dkdT /on 

P'^Pdt ~ ^ \dx^ dy^J dxdx dy dy 

Giving the spacewise as well as time wise temperature gradients as the input to the eq. 

(2.1) will result in the conductivity value ’V as output, ’k’ may or may not be uniform. The 

next article shows the details of the finite-difference scheme. 


2.3 Discretization of GDE 

Equation (2.1) is discretized at each grid point (i,j) and at present time, ’p’. In this 
case, all the space derivatives have been discretized using suitable second order difference 
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schemes. A square grid of Ax = Ay has been used. See Fig 2.2 
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2.3.5 Upper-Left Corner Point (1, n) 
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2.3.6 Interior Grid Points : i=2,n-l and j=:2, n-1 
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2.3.7 Right Boundary : i=n and j=:2, n-1 
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2.3.8 Top Boundary : i=2, n-1 and j=n 


T’P+I T’P 

qC ibA — ibi 
P^p At 


J. i Tr».,-iK,+TU, I -^^.-3+<.-2-5rf„.+2r;a 

''■•n (A)‘ + (SJr / 

2Ax /\ 2Ai J 


+ 2Ai 


(2.9) 


2.3.9 Upper-Right Corner Point (n, n) 
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2A Method of Solution 


Gauss-Seidel iterative method with successive under-relaxation was used to solve the 
set of linear* simultaneous algebraic equations arising out of the discretization of eq.(2.i) at 
each grid point in the computational domain.The solution gives the conductivity value, ’k’ 
at each grid point. The accuracy of the solution depends on grid spacings, time increment 
and the times at which the temperature profiles are taken. The values of p and Cp used 
for obtaining simulated temperatures are 7801 kg/m^ and 473 J/kg K respectively.lt was 
found that a minimum of 16 grid points are required to get a reasonably accurate solution 
for this case since the material ( Mild steel ) used in this study has the conductivity of 45 
W/m.K which was a prior known.The double-precision arithmetic was used. The rate of 
convergence was quite fast The G-S iterative procedure was programmed in Fortran 77. The 
computations were performed using HP-9000 mainframe computer system of I.I.T. Kanpur. 
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Fig.2.1 The Physical problem and the Computational domain 



Fig. 2. 2 The Computational domain showing interior and comer grid points 
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Chapter 3 


Experimental Techniques and Results 


Experimental facilities and techniques used for the transient temperature measure- 
ments in a 4"x4"x0.125" specimen of mild steel to estimate its thermal conductivity are 
first described in this chapter.This is followed by the discussion of the results based on ac- 
tually measured data and after filtering the data.The reasons for the deviation from literature 
value of thermal conductivity of mild steel are given at the end of this chapter. 


3.1 Thermocouples 

A thermocouple is a small electric device used for measuring temperatures. It consists 
of two wires of dissiniilai" metals joined at their ends, forming two junctions. If the junctions 
are at different temperatures, a voltage will be produced that is approximately proportional 
to the temperature difference, a plienomenon discovered by the German physicist Thomas 
Seebeck. 

All temperature measurements in this study were taken with 26 SWG, Type J, insu- 
lated thermocouples, red one indicating iron and white one indicating constantan ( Fig.3 1 
). According to Rangan et al. (1997), Type J thermocouples are made of iron as the positive 
thermocouple element versus constantan ( 51% Cr, 43% Ni ) as the negative thermocou- 
ple element. They are generally recommended for use over the temperature range of -40^C 
to SOO^G. The accuracy is usually ±2% in the range of -40 to 315°C, which the same is 
±0.5%in the range of 375 to 800"C. The Type J thermocouples are inexpensive, mechan- 
ically strong but rapidly deteriorates above 6O0'’C7.In the present experiment temperature 
never exceeded 200'’C. 
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3.2 Formation of Hot junction 


The thennocouple hot junctions were made by flaming in the Glass blowing Workshop, 
I I.T. Kanpur. Prior to that the thermocouple ends were cleaned with the help of emery 
paper.The quality of the junctions were determined by applying certain tension. 


3.3 Caiibration of Thermocouples 

Temperature-measuring devices should be calibrated periodically to evaluate their per- 
formance. Such calibration can be carried out either by comparison with a standard device 
whose accuracy is known. This method has the advantage that the test can be performed at 
any desired temperature range of interest, in a bath. The attainable accuracy depends on the 
uniformity of the temperature distribution in the lest bath and the quality of the reference 
standards used. The standard liqnid-in-glass thermometer can be used as a reference device 
For tests in the range of -100 to 600°C', it is preferable to use an electrically-heated liquid 
bath. Since silicone oil baths are suitable for use as bath liquid in the temperature range 
of 50 to 250‘’C, all the 16 thermocouples used in this experiment were calibrated using an 
electrical heated silicone oil bath with a mercury-in-glass thermometer as the reference de- 
vice. The calibration was done in the Refrigeration and Air-conditioning laboratory of I.I.T 
Kanpur. 


3.4 Heater 

The square specimen of Mild Steel was heated by electric 220/230 Volts AC “Heating 
Mantles” having energy regulator in the range of 0-100 W. The glass fibres enclosed in 
aluminum housing provides high thermal efficiency on the inner surface. The mantles are 
fitted with ON/OFF switches. The maker of the heater is Electro-Mech Instruments, Madras 
See Fig.3.1 for the photograph of the heater. 


3.5 Preparation of the Sample and Assembly 

A mild steel specimen of 4"x4"was cut from a large sheet of 0.125” thickness in the 
Central Workshop, I.I.T.Kanpur. The top and bottom faces were smoothed by polishing and 
the sides were smoothed by grinding. Next, the 16 equally spaced grid points were marked 
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on the top face. At the designated points micro holes were made so that the thermocouple 
hot junctions can be attached using a suitable joining technique. The distance between each 
thermocouple is |inches. In the present case, after experimenting with various joining tech- 
niques such as use of high temperature cement, spot welding etc., it was found that brazing 
gave the best joint. The brazing was carried out in the TA202 Metallurgy laboratory of I.I.T. 
Kanpur. 

After the thermocouples were brazed to tlie designated grid points, the strength of the 
joints were tested by using certain tension, the next step was to connect the cold junctions 
to the connectors by screwing which in turn ai-e attached to two selector switches ( Fig.3 2 
) each capable of handling 10 thermocouples. Finally, the positive and negative ends of the 
last switch were connected to the corresponding ends of the digital temperature recorder. It 
was ensured that the thermocouple number and the number on the switches tally. 

The temperature recorder ( Ftg.3.2 ) has the range of ambient to 200^C. Needless to 
say that the recorder has ambient compensation built into it, The maker of the recorder is 
BlueBell, Kanpur, It has an' LCD display with one digit after the decimal.The accuracy of 
the recorder according to the supplier manual is ±0.5%. This was also checked against a 
mercury-in-glass thermometer reading of ambient as well as boiling water temperature. 


3.6 General Experimental Set-up and Procedure 

The Mild steel square sample supported by four tiny props was placed on an iron screen 
which in turn was positioned on the heater ( Fig. 3. 3 ). The heating mantle was then turned 
on at 50W setting for 28 minutes, 70W setting for 15 minutes and 90W setting for 2 minutes 
till a maximum temperature of 1 86.0 was shown on the digital display of the temperature 
recorder. The overall heating time was 45 minutes. The heating was essentially by radiation 
Table 3.1 at “0” time shows the temperature at the designated grid points at this point in time, 
Cleai'ly, uniform temperature could not be achieved. This is in contrast with the assumption 
in the numerical simulation that the specimen is heated to a uniform temperature. 

The next stage is the cooling of the sample. To do this the hot body was held carefully 
by a pair of tongs and transferred to a cooler surrounding ( at a temperature of 17.2‘’C ). The 
sample was supported by four aluminum props ( with sharp upper ends ) as shown in Fig.3 4. 
During cooling, the temperature at each grid point was recorded every five minutes till 45 
minutes have elapsed. The temperature history during cooling is shown in Table 3.1. The 
entire experiment was carried out in the Fleat Transfer Laboratory. Fig. 3.5 and Fig. 3. 6 show 
the actual thermocouple numbering and corresponding grid point numbering using double 
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subscript notation. 


3.7 Results based on Actually Measured Temperatures 

Table 3.2 shows the conductivity matrices calculated based on the actually measured 
temperatures at successive times starting from time zero. The results clearly are unsatisfac- 
tory as negative conductivity values or not a number (NaN) are not acceptable. However, 
this is not unexpected as it is well-known that inverse problems are highly sensitive to mea- 
surement errors in the input data. Even, 5% measurement errors can produce 50% errors in 
the output. 


3.8 Results based ou Filtered Temperatures 

As mentioned in Art.3.7., the error into the results is usually greater than the error into 
the input data and the solution may be oscillating. Therefore, the solution may be useless 
when the real data are used. In this work a Gram Orthogonal Polynomial method with a 
moving averaging filter window is used for smoothing the noisy data. The method is based 
on a least square appro.vimation. For more details, see Appendix ’B’. Table 3.3 shows the 
transient temperatures after 11 point filtering. Table 3.4 shows the conductivity matrices 
based on the filtered data. Even after filtering, there is virtually no improvement in the results 
pointing to a very high sensitivity of the present numerical algorithm to measurement errors 
The possible resons why even after filtering the results did not improve may be attributed to 
the fact that only timewise filtering was done. Since the present problem is 2D, spacewise 
filtering also should have been performed to get better results. 


3.9 Sources of Experimental Errors 

Since the present algorithm is very sensitive to measurement errors, it is essential that 
the temperature measurements are as accurate as possible. There ai'e several reasons for not 
obtaining the exact temperatures. Some of them are listed below. 

1. In the theoretical formulation, two-dimensional heat conduction is assumed. In actual 
practice, perfect two-dimensionality is very difficult to achieve. 


13 



2. There are 1 6 thermocouples used to measure temperatures at various grid points. These 
sensors are quite long and also thin. So, it produces fm effect thus decreasing the 
accuracy of temperature measurement. 

3. The props on which the plate was heated and cooled were made of materials other than 
mild steel and therefore, actual temperatures were not obtained. 

4. The thermocouples were brazed to the mild steel plate. The filler material for brazing 
is a copper-alloy ( copper, zinc and tin ). The filler material may also have distorted 
the temperature field. 

5. The 16 thermocouple reading should have been taken at the same instant. In the present 
case, there was a time lag since it is not humanly possible to record all the temperatures 
at exactly the same instant. 


3.10 Possible ways to avoid Experimental errors 

To improve the experimental set-up, the following is suggested A PC-based transient 
temperature measurement should be taken with the aid of data acquisition card. Use of non- 
invasive technique such as thermal imaging should be looked into. Better heating and cooling 
arrangements are desirable. Instead of Type J, Type K ( Chromel-Alumel ) thermocouple can 
be used for higher accuracy. High precision temperature recorder should be used. 
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Fig 3.4 The Photograph showing the cooling arrangement 
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Fig.3.5 The location and numbers of the thermocouples 



Fig.3.6 The grid point numbers using double subscript notati 
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i able 3 2 Conductivity Matrices based on Actually Measured Temperat 


NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

1.02 -2.30 54 74 -71.94 

-14.77 -17.07 14.35 30.79 

-32.29 -120.14 -321.57 529.78 

-8.82 -15 34 -90.12 146,20 

14 22 46.91 33.17 34.9! 

-36.78 -57.48 13.86 13.45 

-172.33 -551.15 -151.34 -317.92 

-273.26 -783.14 -107.23 -124.38 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

-1.29 -29.79 17.60 -13.19 

-I5J8 -12.05 14.59 12 49 

-76.16 39.31 -28.77 -93.51 

-.27 53.42 1.63 -9.94 

50,56 -27.33 15 28 2.30 

-27.91 -11.05 12.85 12.70 

-822.64 49 20 -21.31 -66.69 

-131.04 51.11 -3.27 -6.13 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

NaN NaN NaN NaN 

-8,67 -20 01 7.73 1.68 

-11.48 -13.64 8.25 11.18 

11.26 7 26 -12-85 -16.66 

4.99 13.31 -5.54 8.03 



j 158, !6 146.55 138.57 147.50 158.05 14.3 261 139.30 141.431 14435 135.34 148.78 1 138,36 i 143.29 146.64 151.6^ 140.60 

















































































































Table3.^ Conductivity Matrices based on Filtered Temperatures 


16 63 

!4 18 

17 28 

25.28 

-52 86 

-40 78 

-32 68 

-17.97 

-75 51 

-214 82 

92.54 

125 32 

-57 20 

-84.91 

-54.23 

-27.73 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

-1 76 

-24.63 

-7.10 

-5.07 

10.76 

96 20 

14.79 

59.91 

-.94 

5.47 

1.60 

4.88 

864 

-9 49 

39 

-8 41 


-1 35 

-113.57 

-18 19 

-32 54 

497 24 

113441 

-541.42 


NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

-6 76 

-32.10 

-11.46 

-13 14 

-1 43 

23 38 

!8 78 

44.10 

-4 79 

-20 81 

-8.76 

-8.79 

3 26 

18 16 

4 76 

29,13 


-5-21 

5 86 

-12,34 

-9.31 

-65.74 

-85 07 

24 93 

-33.67 

-51.01 

5.63 

-93 43 

11.91 

-96 88 

-145.15 

167.36 

-130,49 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

18.13 

-21 66 

8,84 

9 1 1 

-71.00 

87 35 

-12.58 

44.78 

-1.93 

7.20 

1.78 

-1.71 

14.80 

-1 1.60 

2 33 

-9.89 

-81.92 

-27.19 

14.14 

5.27 

646 55 

224.18 

366 21 

512.72 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

10 03 

20,12 

8.79 

10 75 

-.15 

-151 55 

8 62 

-39.58 

18 90 

-63.40 

11 67 

15.84 

1 07 44 

.349.24 

48 57 

75 95 




Chapter 4 


The Effect of Noise: Simulation Studies 


4.1 Introduction 

The previous chapter has dealt with determination of thermal conductivity of a ho- 
mogeneous plate using actual temperature measurements at discrete points and an inverse 
thermal conduction algorithm. It is evident that the estimated values of thermal conductivity 
are quite sensitive to the accuracy of the input temperature values. The errors at differ- 
ent points and at different times are, in general, not correlated. This can possibly lead to 
sharp spatial and temporal gradients in the temperature field, thereby affecting the calculated 
values of thermal conductivity. This results in severe inconsistencies such as prediction of 
negative or inaccurate conductivity values and spatially varying (in a random manner) con- 
ductivity values even when the test specimen is homogeneous. It is therefore essential to 
understand and quantify the effect of noise on such predictions to make such methods useful 
from a practical point of view. However, in an experimental set-up, it is difficult to control 
the level of noise as an independent parameter. An alternative approach is to simulate such 
problems as outlined in Figure 4.1. Initially, transient temperature data may be generated for 
a given problem (such as a plate of a specified material subjected to a certain temperature 
and then allowed to cool in atmosphere) using the direct heat-conduction formulation. Noise 
may be added to such data using randomly generated numbers within specified ranges. The 
resulting data may then be analysed using inverse calculation methods to predict the thermal 
conductivity of the material. The results can subsequently be compared with the conduc- 
tivity value initially used for the direct problem. If the inverse method being used is based 
on discretization of space, such as a finite difference technique, we generally obtain a range 
of k values for the discrete nodal points. In this case we need to compare the goodness of 
the solution based on some objective measure (parameter) which considers the range of L 
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values obtained. In this work, we propose a number of such criteria and evaluate their rela- 
tive merit. A simulation study also allows us to ascertain the utility of various filters before 
applying inverse methods to noisy data. In the following sections, we present the details of 
the simulation procedure employed and the inferences drawn. 


4,2 Problem Definition: Input Parameters 

We consider a homogeneous plate (4"x 4" x 0.125” ) made of steel. The plate is ini- 
tially at temperature 500^C and is suddenly exposed to ambient at time 1=^0. The properties 
of mild steel are specified as: 


Property 

Units 

Value 

Thermal conductivity, k 

WimK 

45.0 

Density, p 

Kg/m^ 

7801.0 

Specific Heat, Cp 

kJjkgK 

473.0 


We calculate theoretical transient temperature data at 16 grid points. The points are chosen 
corresponding to the locations of thermocouple attachment as described in Chapter 3. 


4.3 Direct Heat- Conduction Formulation 

The direct heat-conduction formulation is employed to determine the transient temper- 
ature data. We assume that there is no temperature variation in the plate-thickness direction 
and accordingly solve the 2D heat-conduction equation; 

dT (d'^T d'^T\ 

dt dy^j 

The boundary conditions and the method of solution using a finite difference formulation are 
desenbed in Appendix ’A’. The step-size in time-increment is choosen in accordance with 
stability requirements ( of the numerical algorithni adopted ) described in Appendix A . 


4.4 Simulation of Noise 


In inverse heat conduction problem there are a number of measured quantifies such 
as temperature, time, senser location, and specimen thickness. We assume that the various 
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qu uilitics uc accurately known except the temperature The temperature measurements ate 
assufiicd to contain the major sources of error or uncertainty. Any known systematic effects 
due to calibration errors, presence of sensors, conduction and convection losses or others 
are assumed to be removed to the extent that the remaining errors may be considered to be 
random. These random enors can then be statistically described. 

The direct problem, as described in Section 4.3 yields simulated temperature data for 
the inverse pioblcm. Noise can then be added to tltis data in terms of random numbers within 
specified ranges. Thus, actual temperature T at a grid point may be given by 

where is the mean (exact) temperature and T is the noise (measurement error). If maxi- 
mum noise level is specified as 5 , we obtain 

T = Tr^ + r.S 

where r is a landom number between -1 and +1. 

The algorithms available for generating random numbers generally yield pseudo ran- 
dom numbers. It is therefore impoitant to consider the statistical description of the numbois 
generated by any such routine. 

One thousand such numbers were generated using NAG routine in the range -1 to +1. 
Their distribution is shown in Fig.4.2 and mean value and Standard deviation are 0.009 & 
0.582 respectively. 


4.5 Results and discussion 

It is of interest to consider the effect of noise in simulated temperatures on inverse cal- 
culation of k value.Since we are using a finite difference algorithm, the solution is obtained 
in terms of spatial distribution of k valucs.The goodness of solution may then be quantified 
on the basis of prediction at a given grid-point, or in a more average sense, in terms of a 
parameter which takes into account the entire field solution. We adopt the latter approach 
and consider various alternative parameters. 


4.5.1 Stability Parameters 

An initial input to the direct problem is the thermal conductivity of the plate material 
For the 4x4 grid being used in the present case. 
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■input — 


^th ^th kth kffi 
^th kih k^h kfji 
kth kill kth kfh 
^th kih kth kth 

where /c^/i is the actual (theoretical) thermal conductivity of the plate material. The coiTe- 
sponding distribution of k values obtained on solving the inverse problem is 

^11 ki2 kis ku 

„ ^21 ^22 ^23 ^24 
^output — 

kzi kz2 k^2 ^34 
^41 ^42 ^43 ^44 

where fcy (i=l“4. j=l-4) represents k values obtained at various gnd point. The koutput 
values will range between k,ain ^>0(3 k„iax> The following parameters are now suggested for 
quantifying the goodness of the results. 



In the ideal case i.e./Cma® = kmin = kth, 21, 22 and Z4 are equal to one and Z3=0. Z1 & Z2 
should be considered together to avoid misleading interpretation under special circumstances 
such as kjyj^ax ~ k^h + A A; and kmm ~ kth ~ AA; Or kmin ~ kth- In practice Z1 alone 

without considering Z2 gave reasonable measure of the goodness of solution.i.e. such special 
case was not encountered in any of the simulations. 
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The parameter Z3 is a more objective measure of the goodness of the solution since 
It considers any deviation in least-square sense. Still Z1,Z2 and Z4 are considered in tins 
work to dcmonsti'ate that physically unacceptable (negative ) values may be obtained if noise 
levels aie high. The use of parameter Z3 fails to bring out this point. 

4.5.2 lYuncation and Roued-off errors 

The use of a numerical algorithm ( both direct and inverse problems ) involves trunca- 
tion and rounding-off errors. Therefore, in general ktfi and koutput will not be identical even 
if no noise were added to the simulated temperatures. A series of typical calculations and the 
icsultingZ values for time steps nl(nli=5, tl=l 16,68 sec) and n2 (n2=35, t2=816.76 sec) are 
shown in Tabic 4.1. It is evident that the values of these parameters deviate from the ideal 
values because of truncation and rounding-off eixors.The deviation is typically 3% 


4.5.3 The effect of noise 

Now we investigate the effect of noise on the stability parameters. Similar to Table 4 I , 
we go tlirough a series of calculations except that now we add random noise to the simulated 
temperatures. Since the temporal gradients change with time, we repeat this simulation for 
various times ( ranging from early to late times ). 

We are interested in the following : 

1 . What is the maximum noise level which still gives us reasonable estimate of k value 
(say, within ±10%of kth) based on inverse calculations. 

2. At what stage of cooling is the numerical scheme most tolerant to noise, i.e., whether 
early_time data is preferred over late„time data. The temporal gradients of temperature 
are much steeper during eaiiy time compared to Iate_times. 

In order to study the maximum permissable noise level, the simulation was run a num- 
bei of times with J as a parameter. The value of 6 ranged between 0 and 1 (insteps of 0.01). 
Each simulation requires temperature information at two time instants (time steps nl and 
n2), which was generated using the direct formulation. Simulations were run for various 
pairs of nl and n2 values to encompass a range of temporal-gradient values including early 

and late times. 

Z values were then calculated for each simulation based on the koutput values obtained 
Typical results on the effect of noise on kontput values and Z parameters as a function of 
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/late time and maxitnuiTi noise level (actual noise at a grid point = random number in 
inge -1 to +I X maximum noise) are shown in Table 4,2.! - 4.2,2, Fig. 4.3. In these 
:s, early lime refers to nl=5 and n2 =35 (i.e. tl =116.68 sec & t2 =816.76 sec) and late 
refers to nl =875 and n2 =905 (i,c. tl =20418.97 sec t2 = 21 1 19.05 scc).Notc that 
rst row of temperature data in these tables corresponds to (5 = 0 , i.e. no noise has been 
1 The following trends emerge from comparision of these results. 

From Fig.4.3 it is clear that the solution is close to theoretical k values up to a critical 
value of 5, and beyond it the solution is unstable often predicting physically unac- 
ceptable (negative) values. Z parameters can therefore be considered as parameters 
measuring the stability of solution This critical value can be measured in terms of 
acceptable deviation in Z value. So, for example, we consider the inverse prediction 
to be acceptable if k is within 10% of k^^- This gives Z ^aceptabie) as 0.9 - 1.1 and 
rptdhif) as 0.1414. Fig 4 4 shows typically the determination of critical value 
of 6 based on Z3 parameter. This value will in general, be different at different times 
during the thermal history. Table 4.3 provides the critical 5 values {Smaz) based on Z1 
paiameter for different times. 

Inverse calculations based on early-time data yield better results than late-time data, 
1 e., the numerical scheme can tolerate larger errors at early times. This is presumably 
related to steeper temporal gradients at early times. 

On calculating (based on simulations) the maximum permissable error in temperature 
at various times during the cooling of the specimen, we find ( Fig,4,5.(b) ) a systematic 
variation of with t. This variation is described well in terms of an exponential 
decay function (Table 4.4). In (2) above, it is supposed that early time data yields better 
results presumably due to steeper gradients. To test this hypothesis, it was decided to 
plot w.r.t. Tlie calculation of with The Values of ^ at various times 
were calculated using noise-free data. Fig.4.6 shows this variation, based both on Z1 
and Z3 parameters. It is found that a linear function (Table 4.5) describes the above 

fKT' 

data quite well, and that has larger value for larger 

The simulations were also run for hypothetical materials with different values of ^ ( 
k =30, 45, 60, 75 ; same p and Cp ). In each case, the results were described well 
by exponential decay function versus t, Fig.4.5) and linear function ((iTnQr versus 
f - Fig.4.7) 

That such a correlation is obtained based on numerical simulation studies, merits fur- 
ittenlion.it is of interest to explain these results in terrrrs of a theoretical frame work 
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In Lhc following section, we attempt to explain these trends in terms of a lumped parameter 
analysis. 


4.5.4 Lumped parameter analysis 

Consider a lumped body, P, suddenly exposed to a surrounding media, S CFig.4.8), at 
a temperature Too- The body P loses the heat to the media S by conduction and the media S 
acts as a heat sink. We assume that the boundary between P and S is such that the equivalent 
heat transfer coefficient in the same as thermal conductivity of P. 

The governing equation for change in temperature of body P is 


pCpV 


dt 


kA 


{Too - T) 
L 


(4.1) 


where L is a characteristic length scale and is of the order ^where V is the volume and A is 
Uie surface area of body P, For an initial temperature To and with both k and Too independent 
of time, the solution of the above equation is 


T = T^-{T^- To) exp (4-2) 

I 

where Tq is the temperature of body P at time t=0. 

In the presence of disturbance, T — Tm-kT" and To = Tom + where subscript m 
refers to exact (mean) temperature and subscript {') represents the disturbance in temperature 
Note that the term T' represents the same physical quantity as 6. A different symbol is used 
here to distinguish between the notations for numerical simulation and theoretical (lumped 
parameter) analysis. 

The above equation used in conjunction with noisy data will give reasonable estimates 
of k only up to a critical noise level above this noise level, the numerical algorithm 
becomes unstable. Therefore T^^^ may be identified as the maximum noise level which is 
acceptable for describing the solution in terms of equation 4.2. 

Therefore we have, 


Tm + T^nax — Too ~ (Too 


(Torn + 


kt \ 
"pC^J 


(4.3) 
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also, Tor exact values of tempeiatures (i.e. if there were no noise). 




CO 


(To, 


- Tom) exp 


kt \ 

pCpiy 


subtracting cq.4.4 from eq.4.3 we obtain 


(4.4) 


-^max '^0~77iax^^P 


kt \ 

pCpiy 


(4.5) 


The above equation implies that maximum permissable error in temperature measurement at 
a given time t which still gives reasonable ( acceptable ) values of k, decays exponentially 
with time.This is in agreement with Fig.4.5. Similar trend is obtained for different kt^ values 
(Fig.4.5). Also, from eq.4.5. 


T' 

-^max 

T 

J- o-mux 


= exp 


kt \ 

pC,Ly 


Therefore a master plot between and kt can be obtained which should be independent 

of kif^ value. Fig.4.9 shows such plot on normalised axes based on die results presented m 
Fig.4,5. 

As far as Fig. 4. 10 is concerned, we can depict the theory behind it in the following 
manner. 

From equation ( 4.1 ), 

— • ^ = T - Too (4-6) 

a at 

where a = and K ~ L 

pCp A 

From equation ( 4.2 ), 


T - Too = - (Too - To) exp (-0t) (4.7) 


where (3 

From cq. 4.6 & eq. 4.7 
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dT d 

~ — (Too - To) exp (4.8) 

Fzom equation ( 4.5 ) ztlso 

T' = r'earp (-/ 5 i) (4.9) 

Dividing equation 4.9 by equation 4.8 , we obtain 


or 


t ' 

max 

t ' 

-**0 

L 2 

dT 

dt 

1 

8 

1 

a 



V 


k )\{To-T^)}' dt 


(4.10) 


The above equation implies that should vary linearly with ^ (since the other 
terms in the above equation are constant and is the maximum permissible noise at 

t=0). This is in good agreement with Fig.4.6-4.7. Again, the data for different k values 
is expected to merge if the results plotted on a normalised x-axis representing 
Although the values tend to come closer, (Fig.4.10), some variation is present presumably 
because all the details of a 2D problem are not captui’ed in a lumped parameter analysis. 


4.6 Effect of Filtering 

In this section we consider the effect of filtering on the random errors that we imposed 
on the simulated temperatures.Here we are using the Gram Orthogonal Polynomial method 
for smoothing the noisy data, which we are getting after the addition of the random errors to 
the temperatures.The Gram Orthogonal Polynomial method has been described in Appendix 
B The main idea is to findout that whether we are getting any improved thermal conductiv- 
ity values at llic same level of noise after applying the filtering Icchniquc.For that we have 
chosen two time instants i.e for early time and for late time. At that time instants we got 
the simulated temperatures ( direct calculation ), and then added the error randomly. We then 
applied the filtering technique to improve the noisy data. The details are tabulated in the 
Table4.6. Based on these results it appears that this filtering technique is not so effective 
with the random errors. 
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Fig.4.8 Lumped body P, exposed to the surrounding media S 



Table 4 ! The effect of truncation and round-off 


errors in inverse problems 
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45 
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45 

45 
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Z1 = I 


22 = 1 


23 = 0 

24: 241 = 1 , 242 = I 


p = 7801kg/m 3 Cp = 473 kJ/kgK 


Direct Simulation 

V 


At t = 1 16,68 sec 


490 90 

493 22 

493 22 

490 90 

493.22 

495.61 

495 61 

493.22 

493 22 

495.61 

495 61 

493.22 

490.90 
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493 22 

490 90 



At t = 

816.76 sec 


456.81 

458 96 
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458 96 

46!. 19 
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Inverse Calculation 

V 
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Tabic 4.2.2 The effect of noise level on Z parameters for laic time steps 
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Table 4.3 Maximum permissible values & gradient for different times based on ZI parameter 


S.No. 


Time Steps 
nl n2 


Gradient range 


Maximum 

Permissible 

Error 


1 

5 

35 

0.048685 

- 

0.049168 

0 . 38 - 

2 

35 

65 

0.045183 


0.045638 

0.335 

3 

65 

95 

0.041938 

- 

0.04236 

0.275 

4 

95 

125 

0.038925 

- 

0.039317 

0.265 

5 

125 

155 

0.036129 

- 

0.036493 

0.29 

6 

155 

185 

0.033534 

- 

0.033872 

0.27 

7 

185 

215 

0.031125 

- 

0.031439 

0.245 

8 

215 

245 

0.02889 

- 

0.029181 

0.195 

9 

245 

275 

0.026814 

- 

0.027085 

0.15 

10 

275 

305 

0,024889 

- 

0.025139 

0.185 

n 

305 

335 

0-023101 

- 

0.023334 

0.17 

12 

335 

365 

0.021442 

- 

0.021658 

0.135 

13 

365 

395 

0.019902 

- 

0.020102 

0.13 

14 

395 

425 

0.018472 

- 

0.018658 

0.125 

15 

425 

455 

0.017145 

- 

0.017318 

0.105 

16 

455 

485 

0.015914 

- 

0.016074 

0.13 

17 

485 

515 

0014771 

- 

0.014919 

0.1 

18 

515 

545 

0.01371 

- 

0.013848 

0.095 

19 

545 

575 

0.012725 

- 

0.012853 

0.07 

20 

575 

605 

0.011811 

- 

0,01193 

0.065 

21 

605 

635 

0.010962 

- 

0.011073 

0.085 

22 

635 

665 

0.010175 

- 

0.010277 

0.068 

23 

665 

695 

0.009444 

- 

0.009539 

0.058 

24 

695 

725 

0-008766 

- 

0-008854 

0.052 

25 

725 

755 

0.008136 

- 

0.008218 

0.042 

26 

755 

785 

0.007552 

- 

0.007628 

0.047 

27 

785 

815 

0.007009 

- 

0.00708 

0.036 

28 

815 

845 

0.006506 

- 

0.006571 

0.034 

29 

845 

875 

0.006038 

- 

0.006099 

0.026 

30 

875 

905 

0.005605 

- 

0.005661 

0.024 


42 



Tabie 4.4 ; The various terms describing the exponential equation 
Equation : y = +Acxp('Bt) 


k 

A 

B 



30 

0,52647 

1.117e-04 

0 

0.0002619 

45 

0.39492 

1.106e-04 

-0.01174 

0.0002699 

60 

0 26526 

9.355e-05 

-0.01618 

3.747e-05 

75 

0 13425 

2.147e-04 

0.05592 

0.0001322 


Table 4.5 : The various terms describing the best fit linear equation 
Equation : y = A+Bx 


k 

A 

B 

Mean 

s.a 

30 

-0.0124 

10.4701 

0.99459 

0.01418 

45 

-0.0172 

7.73875 

0.98733 

0.01567 

60 

1 015C-5 

4.91502 

0.99462 

0,006 

75 

-0.00532 

3.61778 

0.94081 

0.01099 
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Chapter 5 


Conclusions and Scope for Future work 


This work presents a finite-difference based numerical algorithm for the inverse deter- 
mination of thermal conductivity in a 2D square domain using actually measured transient 
temperature data. It is found that for mild steel ( a constant conductivity material ) the algo- 
nlhm is highly sensitive to the measurement errors in the input data, a fact not unexpected 
in inverse problems which are basically ill-posed. In addition to the experimental work, 
a detailed sensitivity analysis is also performed in order to find the maximum permissible 
measurement error ( in temperature ) for such problems. Interestingly, the study reveals that 
the maximum permissible measurement error is greater at early times and decreases expo- 
nentially with time. This means that early time temperature profiles should be used as input 
data. However, even at early times maximum allowable error in the temperature measure- 
ment is quite small which explains why the present experimental input data did not give 
correct conductivity value for mild steel. 

Future effort should be directed towards improving the experimental set-up so that 
high precision temperature measurement is possible. A PC-based experimentation through 
data acquisition card is highly recommended. A non-invasive temperature measurement 
technique such as thermal imaging can also be looked into. A better heating and cooling 
arrangement is necessary. 
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APPENDIX A 


Determination of Simulated Temperature profile 

using Finite-Difference 

A.l Governing Differential Equation 

The governing equation for the constant conductivty case is given by 

A.2 Discretization of GDE 

The details of the discretization of eq.(A.I) at each grid point is described next. See 
Fig.2.1 for the physical problem and computational domain and Fig.2.2 for grid point num- 
bering, The difference schemes of first and second order accuracies are given in Appendix 

A.2.1 Botom Boundary : j=l and i=2, n-1 

This surface is convecting heat to the surroundings.Therefore, 

Af = h[T - r„) 

= I (I’m - 

^ rij-i = Tij+i — 2Aj/| (Tij - Toe) 

Discretizing equation (A.l), we get 



(A.2) 

A.2.2 Top Boundary ; i=2, n-1 and j=n 

This surface is convecting heat to the surroundings. Therefore, 


4 « 



f 


■k%=KT~T^) 


^ (T —T \ 

^ 2Ay “■ k 
^^j+i = ^i,j-i "" 2AyJ (Tij - Too) 


Discretizing equation (A.l), we get 


pC,- 




(A.3) 


A.2.3 Left Boundary : i=l andj=2, n-1 

This surface is convecting heat to the surroundings. Therefore, 

fcf = /»(T - T^) 

= I (Ty - Ito) 

=>■ Ti-ij = Ti+ij - 2Aif (Tjj - Tco) 
Discretizing equation (A.l), we get 


rj-ip+l rj^p 

oC Au iiJ 

P^p At 


'=1 ^ (As/) 


(A.4) 


A.2,4 Right Boundary : i=l and j=2, n-l 

This surface is convecting heat to the surroundings. Therefore, 

-k% = h(T - Too) 

=> Ti+ij = Ti-ij - 2Axj {Tij - Too) 
Discretizing equation (A.l), we get 

pC/^ = k 


(A.5) 
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s 


A.2.5 Interior Grid Points : i=2, n-1 and j =25 n-1 


Discretizing equation (A. I), we get 


rpP+ I 

ZJb2 liii 

P^p At 




(A.6) 


A.2.6 Lower-Left Corner Point : (1,1) 


Discretizing equation (A.l), we get 




'p At 


(A.7) 


A. 2.7 Upper-Left Corner Point : (1, n) 


Discretizing equation (A. I ), we get 

T'p+l T^p ( rpp _OT'P A.'V^ ^OAfhfTP _T ^ TP _9TP XTP OA.i^/^T’P 

■‘i.J -'t.i h J ■‘oo; I 

P^P At Ia^ 

(A.8) 


A.2.8 Lower-Right Corner Point : (n. 1) 


Discretizing equation (A.l), we get 


P^P At \ (A. 




(A.9) 


A.2.9 Upper- Right Corner Point ; (n, n) 


Discretizing equation (A.l), we get 

,r I'ft'-r;, \TU-^Tr^TU,-2Axi{Tf,-T^) r>..-2Tt,-fT^_.-2A,|fe-r~] 
P^p -K < -r 


(A, 10) 
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CHECKING FOR STABILITY 

The coefficient of is taken from each of the equations ( A.2 ) - ( A. 10 ) and is 
it must be greater or equal tozero. In the process, we get diffemet Af's .To ensure 
we choose the minimum of the Af's arising out of the stability condition. 
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Appendix B 

Digital Smoothing Filter : Gram Orthogonal 

Polynomial Method 

In this work a Gram orthogonal polynomial method ( Al-Khalidy, 1998 ) with a moving 
averaging filter windowsis used for smoothing the noisy data. This method is based on a least 
square approximation. This method does not need any information about the beginning and 
the end of the process. This is suitable for use in on-line methods of analysis. Application 
of digital filtering to a series of equally spaced seven data points is shown in fig.B.l. .Digital 
filter replaces each data value T, by a combination of itself and a number of adjacent nodes, 
Thus, 


fiiith) = E bjPj[k) 

where /(ffc) is the smoothed value of the measured temperature at t=k. Let L is the number 
of points used to the left ( past temperatures ) and to the right ( future temperatures ) of the 
central point at time k. The subscript ’n’ refers to the total number of data points used for 
the smoothing process ( n=5 for L=2, n=7 for L=3, etc.). The parameters in eq. (B.l) can be 
calculated as 


- (2Z,+.)4-1)!(2£-j)! 

(2j+l)[(2L)!]" 

(m!)2(24H 

where T denotes the measured values of temperature. Thus, f(t) for each time step as the av- 
erage from Ts^t-L) to T^^t+L) is calculated. This is sometimes called moving window average. 
Notice that — x{x — l){a: — 2) (a: — m-f 1) and m = 1, 2, 3, 

At the first steps the past temperatures can be set to be equal to the initial temperature. 
Similarly, at the last steps the future temperatures can be set equal to the final temperature. 
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The use of 7 to 11 data points is enough to obtain a good approximation. The tempi 
variation between two measurements should be large enough to see variation of tempi 
and should be greater than measurement errors. 
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Appendix C 

Difference Schemes Of First and Second Order 

Accuracy 


C.l Forward difference with error 0{Ax) 


yi ”■ Ax 

^ yi+2-2y^i+j /£ 

{Ax)^ 


C.2 Backward difference with error 0(Aa;) 


7/ — Vi-yj-i 
yz ■“ Ax 

i/£-2yt-i+yi-2 

Vi ” — (Aip 


C,3 Forward difference with error 0{Axy 


J — -.Vi+2+4y,+l-3j/i 

Vi 2 Ax 

<< _ --?/t+3+%4-2-5T/t+i+2yi 

Vi - (Ax)^ 


C.4 Backward difference with error 0(Aa;)^ 


i/i 2 Ax 

^ _ 2?/i-5yf-i+4y,-2-yi'~3 

Vi (Ax)^ 


(c.l) 

(C.2) 


(C.3) 

(c.4) 


(c.5) 

(c.6) 


(C.7) 

(c.8) 


C.5 Central difference with error (9(Aa:)^ 


